%-------------------------------------------------------------------------;
% Horizon Bias and the Term Structure of Equity Returns;
% This code replicates:
%   - Tables 1, 2, 3, and 6
%   - All figures
%-------------------------------------------------------------------------;

clear all;
clc;
close all;

%-----------------------------;
% Import term structure data;
%-----------------------------;

filename='Term_structure_variables';
data_ts=xlsread(filename);
b=size(data_ts,1); 

Div_ret=data_ts(2:end,6);
SP_ret=data_ts(2:end,4);

Div_ret=log(1+Div_ret);
SP_ret=log(1+SP_ret);
                                                 
etp=SP_ret-Div_ret;

SP_price=data_ts(1:end,1);
SP_Div=data_ts(1:end,3);
SP_strip_18=data_ts(1:end,5);

SP_pd=log(SP_price./SP_Div);
SP_pd_strip=log(SP_strip_18./SP_Div);

SP_pd=SP_pd(1:b-12);
SP_pd_strip=SP_pd_strip(1:b-12);

%-------------------------------;
% Import explanatory variables;
%-------------------------------;

filename='Explanatory_variables';
data_ev=xlsread(filename);

data_ev=data_ev(1:b-12,:);

year=data_ev(1:end,1);
month=data_ev(1:end,2);

Ltg=data_ev(1:end,3);
Stg=data_ev(1:end,4); 

%Take logs;
Ltg=log(1+Ltg);
Stg=log(1+Stg);

%Unaffiliated;
Ltg_un=data_ev(1:end,5);
Stg_un=data_ev(1:end,6); 

Ltg_un=log(1+Ltg_un);
Stg_un=log(1+Stg_un);

%Skewness;
skew1=data_ev(1:end,7);
skew1=(skew1-mean(skew1))./std(skew1);

skew5=data_ev(1:end,8);
skew5=(skew5-mean(skew5))./std(skew5);
Diff_skew=skew5-skew1;

%Uncertainty;
uncert1=data_ev(1:end,9);
uncert1=(uncert1-mean(uncert1))./std(uncert1);

uncert5=data_ev(1:end,10);
uncert5=(uncert5-mean(uncert5))./std(uncert5);
Diff_uncert=uncert5-uncert1;

%Other variables capturing business cycle;
rec=data_ev(1:end,11);
default=data_ev(1:end,12);
term=data_ev(1:end,13);
cay=data_ev(1:end,14);

%Horizon bias;
acorr=autocorr(Stg);
theta_0=mean(Stg);
theta_1=acorr(13,1);

Stg_impl_one=(1/5)*(4*theta_0+3*theta_0*theta_1+2*theta_0*(theta_1^2)+theta_0*(theta_1^3));
Stg_impl_two=(1/5)*(1+theta_1+theta_1^2+theta_1^3+theta_1^4)*Stg;

Stg_impl=Stg_impl_one*ones(size(Stg_impl_one,1),1)+Stg_impl_two;

HB=Ltg-Stg_impl;

%Horizon bias - unaffiliated;

acorr=autocorr(Stg_un);
theta_0=mean(Stg_un);
theta_1=acorr(13,1);

Stg_impl_one=(1/5)*(4*theta_0+3*theta_0*theta_1+2*theta_0*(theta_1^2)+theta_0*(theta_1^3));
Stg_impl_two=(1/5)*(1+theta_1+theta_1^2+theta_1^3+theta_1^4)*Stg_un;

Stg_impl_un=Stg_impl_one*ones(size(Stg_impl_one,1),1)+Stg_impl_two;

HB_un=Ltg_un-Stg_impl_un;

%------------------------------;
%Summary statistics: Table 1;
%------------------------------;

auto_SP_ret=autocorr(SP_ret);
auto_Div_ret=autocorr(Div_ret);
auto_etp=autocorr(etp);
auto_Ltg=autocorr(Ltg);
auto_Stg=autocorr(Stg);
auto_Stg_impl=autocorr(Stg_impl);
auto_HB=autocorr(HB);

auto_SP_pd=autocorr(SP_pd);
auto_SP_pd_strip=autocorr(SP_pd_strip);

auto_term=autocorr(term);
auto_default=autocorr(default);
auto_cay=autocorr(cay);
auto_rec=autocorr(rec);

auto_diff_skew=autocorr(Diff_skew);
auto_diff_uncert=autocorr(Diff_uncert);

Table1=[length(SP_ret),mean(SP_ret),std(SP_ret),min(SP_ret),max(SP_ret),auto_SP_ret(2,1),auto_SP_ret(13,1);
    length(Div_ret),mean(Div_ret),std(Div_ret),min(Div_ret),max(Div_ret),auto_Div_ret(2,1),auto_Div_ret(13,1);
    length(etp),mean(etp),std(etp),min(etp),max(etp),auto_etp(2,1),auto_etp(13,1);
  
    length(Ltg),mean(Ltg),std(Ltg),min(Ltg),max(Ltg),auto_Ltg(2,1),auto_Ltg(13,1);
    length(Stg),mean(Stg),std(Stg),min(Stg),max(Stg),auto_Stg(2,1),auto_Stg(13,1);
    length(Stg_impl),mean(Stg_impl),std(Stg_impl),min(Stg_impl),max(Stg_impl),auto_Stg(2,1),auto_Stg(13,1);
    length(HB),mean(HB),std(HB),min(HB),max(HB),auto_HB(2,1),auto_HB(13,1);
        
    length(Diff_skew),mean(Diff_skew),std(Diff_skew),min(Diff_skew),max(Diff_skew),auto_diff_skew(2,1),auto_diff_skew(13,1);
    length(Diff_uncert),mean(Diff_uncert),std(Diff_uncert),min(Diff_uncert),max(Diff_uncert),auto_diff_uncert(2,1),auto_diff_uncert(13,1);
        
    length(SP_pd),mean(SP_pd),std(SP_pd),min(SP_pd),max(SP_pd),auto_SP_pd(2,1),auto_SP_pd(13,1);
    length(SP_pd_strip),mean(SP_pd_strip),std(SP_pd_strip),min(SP_pd_strip),max(SP_pd_strip),auto_SP_pd_strip(2,1),auto_SP_pd_strip(13,1);
    length(term),mean(term),std(term),min(term),max(term),auto_term(2,1),auto_term(13,1);
    length(default),mean(default),std(default),min(default),max(default),auto_default(2,1),auto_default(13,1);
    length(cay),mean(cay),std(cay),min(cay),max(cay),auto_cay(2,1),auto_cay(13,1);
    length(rec),mean(rec),std(rec),min(rec),max(rec),auto_rec(2,1),auto_rec(13,1);]

%-----------------------;
% Table 2;
%-----------------------;

a=[SP_ret(1:end-12),Div_ret(1:end-12),etp(1:end-12),Ltg(2:end),Stg(2:end),HB(2:end),Diff_skew(2:end),Diff_uncert(2:end)];

A=corr(a);
Table2=A(1:8,3:8)

%Check correlation between long-term growth and differences in prices of long- and short-term assets; 
corr(Ltg,SP_pd-SP_pd_strip);
corr(HB,SP_pd-SP_pd_strip);

%------------------------------;
%Figures;
%------------------------------;

h=12;

%Sum returns over h periods;
ret_SP=filter(ones(h,1),1,SP_ret);
ret_Div=filter(ones(h,1),1,Div_ret);

ret_SP=ret_SP(h:end);
ret_Div=ret_Div(h:end);

x=linspace(-1,1);

%SP 500 return;
figure(1)
plot(ret_SP,'-','Color',[0.15,0.15,0.7],'LineWidth',1)
hold on;
plot(zeros(size(ret_SP,1),1),'--','Color',[0.25,0.25,0.25],'LineWidth',1)
legend('S&P 500 return')
set(gca, 'XTick', 6:24:length(year));
set(gca, 'XTickLabel', year(6:24:length(year)));
axis([1 length(ret_SP) -1.1 1]);
set(gca,'FontSize',18)

%Dividend strip return;
figure(2)
plot(ret_Div,'-','Color',[0.15,0.7,0.15],'LineWidth',1)
hold on;
plot(zeros(size(ret_SP,1),1),'--','Color',[0.25,0.25,0.25],'LineWidth',1)
legend('Dividend strip return')
set(gca, 'XTick', 6:24:length(year));
set(gca, 'XTickLabel', year(6:24:length(year)));
axis([1 length(ret_SP) -1.1 1]);
set(gca,'FontSize',18)

%Realized equity term premium;
figure(3)
plot((ret_SP-ret_Div),'-','Color',[0.7,0.15,0.15],'LineWidth',1)
hold on;
plot(zeros(size(ret_SP,1),1),'--','Color',[0.25,0.25,0.25],'LineWidth',1)
legend('Realized equity term premium')
set(gca, 'XTick', 6:24:length(year));
set(gca, 'XTickLabel', year(6:24:length(year)));
axis([1 length(ret_SP) -1.1 1]);
set(gca,'FontSize',18)

%Long-term growth;
figure(4)
plot(Ltg,'-','Color',[0.15,0.15,0.7],'LineWidth',1)
hold on;
plot(zeros(size(ret_SP,1),1),'--','Color',[0.25,0.25,0.25],'LineWidth',1)
legend('Long-term growth rate')
set(gca, 'XTick', 6:24:length(year));
set(gca, 'XTickLabel', year(6:24:length(year)));
axis([1 length(ret_SP) 0.05 0.2]);
set(gca,'FontSize',18)

%Short-term growth;
figure(5)
plot([Stg],'-','Color',[0.15,0.7,0.15],'LineWidth',1)
hold on;
plot(zeros(size(ret_SP,1),1),'--','Color',[0.25,0.25,0.25],'LineWidth',1)
legend('Short-term growth rate')
set(gca, 'XTick', 6:24:length(year));
set(gca, 'XTickLabel', year(6:24:length(year)));
axis([1 length(ret_SP) -0.31 0.31]);
set(gca,'FontSize',18)

%Horizon bias;
figure(6)
plot(HB,'-','Color',[0.7,0.15,0.15],'LineWidth',1)
hold on;
plot(zeros(size(ret_SP,1),1),'--','Color',[0.25,0.25,0.25],'LineWidth',1)
legend('Horizon bias')
set(gca, 'XTick', 6:24:length(year));
set(gca, 'XTickLabel', year(6:24:length(year)));
axis([1 length(ret_SP) -0.05 0.08]);
set(gca,'FontSize',18)

%Horizon-differenced skewness;
figure(7)
plot([skew5-skew1],'-','Color',[0.7,0.15,0.15],'LineWidth',1)
hold on;
plot(zeros(size(ret_SP,1),1),'--','Color',[0.25,0.25,0.25],'LineWidth',1)
legend('Horizon-differenced skewness')
set(gca, 'XTick', 6:24:length(year));
set(gca, 'XTickLabel', year(6:24:length(year)));
axis([1 length(ret_SP) -2 2]);
set(gca,'FontSize',18)

%Horizon-differenced uncertainty;
figure(8)
plot([uncert5-uncert1],'-','Color',[0.7,0.15,0.15],'LineWidth',1)
hold on;
plot(zeros(size(ret_SP,1),1),'--','Color',[0.25,0.25,0.25],'LineWidth',1)
legend('Horizon-differenced uncertainty')
set(gca, 'XTick', 6:24:length(year));
set(gca, 'XTickLabel', year(6:24:length(year)));
axis([1 length(ret_SP) -2 2]);
set(gca,'FontSize',18)

%Horizon bias and realized equity term premium;

%Set dates such that year will fit middle of each year;
dates=(1996:(1/12):2016.999)';
dates=dates-(5/12);

figure(9)
[AX,H1,H2] = plotyy(dates,HB,dates,(ret_SP-ret_Div));
ylabel(AX(1),'Horizon bias','FontSize',14) % left y-axis;
ylabel(AX(2),'Realized equity term premium','FontSize',12) % right y-axis;
set(AX,{'ycolor'},{'b';'[0,0.5,0]'})  % Left blue, right dark green;
set(AX(1),'YLim',[-0.06 0.09])
set(AX(1),'YTick',(-0.06:0.03:0.09))
set(AX(2),'YLim',[-1.2 0.8])
set(AX(2),'YTick',(-1.2:0.4:0.8))
set(AX(1),'xLim',[1995.5 2016.5])
set(AX(2),'xLim',[1995.5 2016.5])
set(AX(1),'xTick',(1996:2:2017))
set(AX(2),'xTick',(1996:2:2017))
set(AX(1),'FontSize',18)
set(AX(2),'FontSize',18)
set(H1,'Linewidth',1,'Color','b');
set(H2,'Linewidth',1,'LineStyle','--','Color','[0,0.5,0]');

%--------------------;
%Table 3;
%--------------------;

y=ret_SP-ret_Div;
a=HB;
f=size(y,1);

b=find(a(:)>=mean(a));
c=find(a(:)<mean(a));

d=y(b,:);
e=y(c,:);

Table3(1:2,1)=[mean(d);mean(e)];

%Unaffiliated;

a=HB_un;
f=size(y,1);

b=find(a(:)>=mean(a));
c=find(a(:)<mean(a));

d=y(b,:);
e=y(c,:);

Table3(1:2,2)=[mean(d);mean(e)]

%--------------------;
%Table 6;
%--------------------;

%Skewness;

a=HB.*(skew5-skew1);
f=size(y,1);

b=find(a(:)>=mean(a));
c=find(a(:)<mean(a));

d=y(b,:);
e=y(c,:);

Table6(1:2,1)=[mean(d);mean(e)];

%Uncertainty;

a=HB.*(uncert5-uncert1);
f=size(y,1);

b=find(a(:)>=mean(a));
c=find(a(:)<mean(a));

d=y(b,:);
e=y(c,:);

Table6(1:2,2)=[mean(d);mean(e)]

%-----------------------;
% Figure 5;
%-----------------------;

a=HB;
f=size(y,1);

for prc=1:101;
b=find(a(:)>=prctile(a,prc-1));
c=find(a(:)<=prctile(a,prc-1));

d=y(b,:);
e=y(c,:);

Table_hb(prc,1:2)=[mean(d) mean(e)];
end;

Table_hb=Table_hb(1:10:91,:);

y_hb=[Table_hb(:,1)];

x=[0:10:90];
constant=ones(10,1).*(mean(y));

figure(10)
bar(x,[y_hb]);
legend('Horizon bias');
ylabel('Mean ETP')
xlabel('Percentile')
ylim([-0.5 0.1]);
xlim([-5 95]);
set(gca,'FontSize',18)

%---------------------------------------------------------;
% Figure 6: Interactions with skweness and uncertainty;
%---------------------------------------------------------;

%HB*(Long-term minus short-term skewness);

a=HB.*(skew5-skew1);
f=size(y,1);

for prc=1:101;
b=find(a(:)>=prctile(a,prc-1));
c=find(a(:)<=prctile(a,prc-1));

d=y(b,:);
e=y(c,:);

Table_hb_skew(prc,1:2)=[mean(d) mean(e)];
end;

Table_hb_skew=Table_hb_skew(1:10:91,:);
y_hb_skew=[Table_hb_skew(:,1)];

%HB*(Long-term minus short-term uncertainty);

a=HB.*(uncert5-uncert1);
f=size(y,1);

for prc=1:101;
b=find(a(:)>=prctile(a,prc-1));
c=find(a(:)<=prctile(a,prc-1));

d=y(b,:);
e=y(c,:);

Table_hb_uncert(prc,1:2)=[mean(d) mean(e)];
end;

Table_hb_uncert=Table_hb_uncert(1:10:91,:);
y_hb_uncert=[Table_hb_uncert(:,1)];

x=[0:10:90];
constant=ones(10,1).*(mean(y));

figure(11)
bar(x,[y_hb_uncert y_hb_skew]);
legend('Horizon bias*(Horizon-differenced uncertainty)','Horizon bias*(Horizon-differenced skewness)');
ylabel('Mean ETP')
xlabel('Percentile')
ylim([-0.5 0.20]);
xlim([-5 95]);
set(gca,'FontSize',18)

%---------------------------------------;
%Bias ve. rational component;
%---------------------------------------;

Ltg=data_ev(1:end,15);
Ltg_r=data_ev(1:end,16);
Ltg_bias=data_ev(1:end,17);

Ltg=log(1+Ltg);
Ltg_r=log(1+Ltg_r);
Ltg_bias=log(1+Ltg_bias);

Stg=data_ev(1:end,18);
Stg_r=data_ev(1:end,19);
Stg_bias=data_ev(1:end,20);

Stg=log(1+Stg);
Stg_r=log(1+Stg_r);
Stg_bias=log(1+Stg_bias);

%Rational component of horizon-differenced earnings;

a=autocorr(Stg_r);
theta_0=mean(Stg_r);
theta_1=a(13,1);

Stg_impl_one=(1/5)*(4*theta_0+3*theta_0*theta_1+2*theta_0*(theta_1^2)+theta_0*(theta_1^3));
Stg_impl_two=(1/5)*(1+theta_1+theta_1^2+theta_1^3+theta_1^4)*Stg_r;

Stg_impl_r=(Stg_impl_one*ones(size(Stg_impl_one,1),1)+Stg_impl_two);

HB_r=Ltg_r-Stg_impl_r;

%Biased component of horizon-differenced earnings;

a=autocorr(Stg_bias);
theta_0=mean(Stg_bias);
theta_1=a(13,1);

Stg_impl_one=(1/5)*(4*theta_0+3*theta_0*theta_1+2*theta_0*(theta_1^2)+theta_0*(theta_1^3));
Stg_impl_two=(1/5)*(1+theta_1+theta_1^2+theta_1^3+theta_1^4)*Stg_bias;

Stg_impl_bias=(Stg_impl_one*ones(size(Stg_impl_one,1),1)+Stg_impl_two);

HB_bias=Ltg_bias-Stg_impl_bias;

%--------------------------;
%Figure 7;
%--------------------------;

%Biased component;

a=HB_bias;

f=size(y,1);

for prc=1:101;
b=find(a(:)>=prctile(a,prc-1));
c=find(a(:)<=prctile(a,prc-1));

d=y(b,:);
e=y(c,:);

Table_ltg(prc,1:2)=[mean(d) mean(e)];
end;

Table_ltg=Table_ltg(1:10:91,:);

y_hb_bias=[Table_ltg(:,1)];

x=[0:10:90];
constant=ones(10,1).*(mean(y));

figure(12)
bar(x,[y_hb_bias]);
legend('Biased component');
ylabel('Mean ETP')
xlabel('Percentile')
ylim([-0.4 0.1]);
xlim([-5 95]);
set(gca,'FontSize',18)

%Rational component;

a=HB_r;
f=size(y,1);

for prc=1:101;
b=find(a(:)>=prctile(a,prc-1));
c=find(a(:)<=prctile(a,prc-1));

d=y(b,:);
e=y(c,:);

Table_ltg(prc,1:2)=[mean(d) mean(e)];
end;

Table_ltg=Table_ltg(1:10:91,:);
y_hb_rational=[Table_ltg(:,1)];

x=[0:10:90];
constant=ones(10,1).*(mean(y));

figure(13)
bar(x,[y_hb_rational]);
legend('Rational component');
ylabel('Mean ETP')
xlabel('Percentile')
ylim([-0.4 0.1]);
xlim([-5 95]);
set(gca,'FontSize',18)





